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Abstract Simultaneous deterministic dynamics of multiple populations de- 
scribed by a large set of ODE’s is considered in the phase space of population 
sizes and ODE’s parameters. The problem is formulated as a multidimensional 
phase-space conservation law and is solved explicitly for non-interacting multi¬ 
population models. Solutions for populations competing for a limited resource, 
populations with migration, as well as populations with phase-space interac¬ 
tions are obtained by simple iterative methods. 
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1 Introduction 

Certain natural objects and phenomena, not immediately associated with pop¬ 
ulation dynamics, can still be described as a large number of simultaneously 
developing populations. For example, forests consist of plants that in turn 
consist of many cells. Also cells themselves contain significant populations of 
mitochondria - small organelles capable of multiplying on their own. Thus, 
considering the mitochondria of a single cell to be the elementary population, 
one can view the cells of a plant as a population of populations, or a multi¬ 
population. Alternatively, viewing a single plant as an elementary population 
of cells, a field of plants is a population of populations, etc. (see, e.g., Merchant 
et al. 1960; Carrie et al. 2012). 

Although, classical population dynamics is a very well developed subject 
(see, e.g., Shonwinkler and Herod 2009; Begon et al. 2006; Newman et al. 
2014), its main concern has traditionally been the growth of a single or at most 
of a few populations simultaneously. The question of dynamics of dispersing 
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populations has led to the development of the metapopulation theory (e.g., 
Levin 1969; Eriksson et al. 2014) that treats populations with spatially disjoint 
habitats as a single metapopulation, and reduces the mathematical description 
to a single ODE. Alternatively, it is sometimes advantageous to consider the 
space-time dynamics of all members of all populations simultaneously, perhaps, 
subdividing them in a few species as well. This leads to spatially-resolved PDE 
models, such as the Fischer equation (Fischer 1930), that have produced many 
valuable insights over the years (see, e.g., Edelstein-Keshet 2005). 

Studying multiple populations with large systems of ODE’s or spatially- 
resolved PDE’s does not always provide explicit answers to the pertaining ques¬ 
tions. For example, in agriculture and ecology (e.g., Hautala and Hakojarvi 2011; 
Lin et al. 2013; Xu et al. 2014; Velazquez et al. 2014), where the growth of 
an individual plant or species is typically described by a logistic or a reaction- 
type ODE, the actual question of interest or the measured data is often the 
histogram of plant or population sizes. Predicting the time-evolution of such 
a histogram is not a trivial problem, since not only the initial conditions but 
also the various parameters of the ODE may differ from plant to plant. An¬ 
other complication arises if populations interact in some way, e.g., there may 
be a competition for a limited resource or migration of individuals from one 
population to another. Populations may also interact on a more esoteric level, 
e.g., their growth rates may depend on the distribution over the phase space. 

Histograms or their continuous analogue - phase-space distributions - are 
typically associated with random variables and stochastic processes and are 
less common in deterministic setting. Nonetheless, the temporal evolution of 
such deterministic distribution functions may be useful in the analysis of bio¬ 
logical data, since it may help to distinguish deterministic causal effects from 
stochastic phenomena and random correlations. There is nothing new in the 
phase-space analysis of deterministic dynamics as such. In fact, it is often 
performed on the level of single trajectories (phase-space curves), e.g., the 
predator-prey system (see, e.g., Edelstein-Keshet 2005). The novelty of the 
present contribution is in the systematic extension of these ideas from single 
trajectories to distributions over the phase space. 

To predict the time evolution of the multi-population distribution function 
we employ an n-dimensional phase-space conservation law, where the mathe¬ 
matical form of the phase-space current is determined by the dynamic equa¬ 
tions (coupled or uncoupled ODE’s) of individual populations. This approach 
has been inspired by the Lifshitz-Slyozov-Wagner method from materials sci¬ 
ence (Lifshitz and Slyozov 1961; Wagner 1961; Kampman and Wagner 1991; 
Myhr and Grong 2000; Collet and Goudon 2000; Laurencot 2002). Concep¬ 
tually it may be viewed as a version of the Fokker-Planck equation tuned, in 
our case, to the specific needs of the multi-population dynamics. The present 
approach is also similar in spirit to the Liouville equation of Hamiltonian me¬ 
chanics, except that the population dynamics is not necessarily Hamiltonian. 
The closest analogue in the field of population dynamics is, probably, the 
state-space method discussed in (Newman et al. 2014). 
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To illustrate the benefits and limitations of the proposed framework we in¬ 
vestigate a range of population models relevant to plant biology and migration 
studies. We focus on problems that have either completely explicit solutions or 
can be solved with a simple iterative algorithm. In the subsequent sections we 
consider populations with unlimited (exponential) growth, populations whose 
growth is limited by the carrying capacity of the medium (logistic equation), 
populations competing for a single limited resource (e.g. oxygen, food, etc.), 
populations with migrating members, and populations whose rate of growth 
depends on the distribution function. Whenever possible we include both the 
population size and the parameters of ODE’s in a single multidimensional 
distribution evolving in time. 


2 Populations with unlimited growth 


Let a very large number of populations be growing in accordance with the 
following equations: 

din ■ 

—^=airii, m( 0) > 0, i = (1) 

dt 

where oti > 0 are the growth constants. The variable rij will be referred to as 
the size of the *-th population, but, in practice, it may denote many different 
things. For example, it can be the length, volume, or the biomass of a plant. 
It can also be a continuous approximation of the number of animals within a 
habitat, or the number of cells in a plant, or the number of mitochondria in a 
cell. The units of time are also completely arbitrary here. 

Instead of solving all these ODE’s we consider the time evolution of the 
distribution function u(n, a , t), whose integral over the size variable n and the 
growth constant variable a represents the total number of populations, which 
is conserved: 


OO OO 

J J u(n, a, t ) dn da = P. 


o o 


( 2 ) 


Integrating the distribution over a only, 


OO 

J u(n,a,t) da = p(n,t), (3) 

o 

one obtains a continuous approximation of the distribution of populations over 
their sizes at time t. The total size N(t) of all populations at time t can be 
obtained as: 
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We assume that the sufficiently smooth continuous approximation of the 
initial distribution u(n,a, 0) = uo(n,a) is given. The latter is interpreted 
as a two-dimensional distribution of populations over their initial sizes n o 
and their growth constants a. Since the dynamics of all populations is purely 
deterministic, the functions uq and u are not probability density functions. 
However, upon suitable normalization, they can be interpreted as such for a 
randomly selected population. 

The dynamics of populations causes the distribution function u to change 
in time. To account for these changes we introduce a phase-space current with 
the density J. As a consequence of the conservation law © the distribution 
function u obeys the phase-space continuity equation: 

| + V.J = 0, (5) 

where the vector J, in general, has as many components as there are phase- 
space variables, i.e. two in the present case: 


3 = (uv n ,uv a ). (6) 

The corresponding ‘velocities’ are deduced from the dynamic equation (JT|) , 
i.e., in the present case: v n = n' = an and v a = cd = 0. Since J has only 
one nonzero component, the continuity equation ([5| reduces to the following 
Cauchy problem: 

du du , 

—h an— —I- au = 0, u[n, a, 0) = uq(ti, a). (7) 

at on 

The solution of this problem, obtained with the method of characteristics, is 

u(n,a,t) = e~ at u 0 (ne~ at ,a). (8) 

An example of the time evolution of u(n, a, t) is shown in Fig. [TJ where the ini¬ 
tial distribution u^{n, a) is a Gaussian centered at n c = 15 and a c = 20, which 
was set to zero for a,n < 0 (see the top-left image). The size distributions of 
populations p(n,t) (bottom row of Fig. [l]) were obtained by numerical integra¬ 
tion over a and demonstrate a progressively fattening tail due to the rapidly 
growing populations corresponding to the upper part of u(n,a,t), i.e., larger 
growth rate constants a. An interesting feature in the evolution of u(n, a, t) 
that becomes more pronounced at larger i’s, is the gradual downwards shift 
of the maximum of u(n, a, t) with time. This is due to the fact that the pop¬ 
ulations with smaller a’s tend to grow more evenly (synchronously) even if 
their initial n’s are different, whereas, populations with larger a’s and slightly 
different initial n’s rapidly grow apart. 
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t=0.0 t=0.015 t=0.03 






Fig. 1 Time evolution of u(n, a , t) according to (|8|. Top row - snapshots at t = 0.0, 0.015, 
and 0.03; Bottom row - corresponding numerically integrated p(n, t ) showing the distribu¬ 
tion of populations over their sizes (the horizontal, n-axis, of these plots is shared with the 
images above). 


3 Populations with logistic growth 

Let the dynamics of populations be governed by the logistic equations 
drii 


dt 


= lirii{ki — rii), rii(0) > 0, i = 1,2,..., 


( 9 ) 


where 7 * > 0 is a rate constant, and hi > 0 is a constant determining the max¬ 
imum sustainable population. The distribution function u(n,j,k,t) over the 
three-dimensional (n, 7 , k) phase space satisfies the following Cauchy problem: 

du du 

— +'yn(k - n)— + 7 (fc - 2n)u = 0 , u(n, 7 , k, 0 ) = u 0 (n, 7 , k). ( 10 ) 

Here too the solution can be obtained with the method of characteristics. Since 
it is more involved than in the previous case, we shall provide the intermediate 
steps. The characteristic equations of (10) are 


dt dn du . . 

— = 1 , — = 7 n(fc — n), — = —7 (ft — 2n)u. 

as as as 


( 11 ) 


This leads to the following relations: 


C n = 


| k- 


„-i kt 


C 1 


i\k — 


u = 


( 12 ) 
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t=0.0 t=0.001 t=0.002 



Fig. 2 Time evolution of u(n, 7 , k , t) according to \l6\ with fixed k = 80. Top row - 
snapshots at t = 0.0,0.001, and 0.002; Bottom row - corresponding numerically integrated 
p(n, t) showing the distribution of populations over their sizes. 


We shall restrict ourselves to the segment 0 < n < k, so that \k — n\ = k — n. 
The general form of the solution to (10) in this case is 


f{p) 

n(k — n ) ’ 


P = 


o~l kt 


k — n 


(13) 


The function f(jp ) should be chosen in such a way that the initial condition 
is satisfied. To this end we introduce two smooth functions of p that at t = 0 
behave as follows: 


kp 


p+ 1 


= n, 


k 2 p 


t—O 


= n(k — n). 


t—0 


{p+ l ) 2 

The initial condition from (10) is satisfied if these functions are used as 
k 2 p 


(14) 


n{k — n)(p + l ) 2 


u o 


kp 

p + 1 




= «o(n,7)k)- 


t—0 


Thus, the solution of the Cauchy problem (10) on the strip 0 < n < k is 

fc 2 e -7fei 


i(n,7,fc,i) = 


(■ne lkt + k — n) z 


uq 


kne lkt 

ne -lkt -)- — 


-,7 


(15) 


(16) 


and for n = k, we obviously have u(k, 7 , fc, t) = uo(k , 7 , k). 

An example of the time evolution of u(n, 7 , k, t) for a fixed k = 80 is 
shown in Fig. |2j where the initial distribution (n-, 7,80) (top left) is the 
same Gaussian as in the example of Fig. [T] As t grows the corresponding 
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size distributions p(n, t) (bottom row) tend to a distribution localized at the 
right boundary n = k. This time the maximum of the distribution u(n, 7 , k, t) 
shifts upwards (eventually), meaning that for larger t's the populations with 
larger growth rates 7 will dominate, as they will rapidly accumulate around 
the maximum attainable size. 


4 Populations competing for a resource 


Consider populations competing for a limited resource c that controls their 
growth rate. This control substance is absorbed equally well by all members of 
all populations. In response, the populations grow. However, we assume that 
the growth of a population slows down and eventually stops if the supply of c 
drops to some minimal value c m ; n , which is here taken to be zero for simplicity. 
The following coupled system of equations describes this situation: 


^ = - 7 X! m, c(0) = c 0 ,7 > 0, 
d n ■ 

“TT = fiiTliC, 7lj(0), Pi > 0, * = 1,2, ... 


(17) 


Integrating the first equation and substituting the result in the second equation 
we arrive at: 


drii 

dt 



i= 1,2,... 


(18) 


The sum inside these equations is the total size N(t) of all populations at time 
t. which in terms of the distribution function u(n, /3, t) can be written as 


OO OO 


^n,(f) = N(t) = 


13, t)dnd0. 


0 0 


(19) 


Thus, the phase-space dynamics of the distribution function will be governed 
by a weakly nonlinear equation: 

du du 

— +/3cn—+pcu = 0, u(n,/3,0) = u 0 (n,/3), (20) 

where 


t 

c(t) = c 0 — 7 J N(t')dt r . ( 21 ) 

0 

The questions of existence and uniqueness of the solutions for this and similar 
equations appearing below are treated in (Collet and Goundon 2000; Lauren- 
cot 2002 ). 
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As expected, the equation (201 confirms that the distribution u stops 


changing in time (becomes stationary) when the supply c is exhausted. Al¬ 
though, this behaviour is asymptotic for t —¥ oo, one may expect that the 
changes in u are negligible for all t > T for some sufficiently large T. 


Let us focus on finding the solution u of the linear Cauchy problem (20) 
for a given function c(t). Exploiting the direct analogy with 
of (20) is found to be 


the solution 


where 


u(n,/3,t)=e P^u,o(ne ^^,/3), 


€{t)= f c(t) dr. 
Jo 


( 22 ) 


(23) 


Equation fl22| is an implicit solution of the original problem, since £(t) depends 
on u via (19), plj ), and (23). 

Let u k (n,/3) = u(n,/3,t k ), N k = N(t k ), c k = c(t k ), £ fc = £(4), and t k+1 = 
t k + At. For sufficiently small At we can approximate the solution with the 
following time-stepping scheme: 

Algorithm 4 

— Given: £ 0 = 0, Co > 0, Uo(n,(3), and a sufficiently large Q. 

— For k = 0,1, 2,..., while c k > 0, do: 

Cfc+i = C/c A C k At, 

u k+1 (n,/3) = e _/35fc + 1 -u 0 (ne _,9?fc + 1 ,/?), 

N k+1 = jj nu k+ i(n,/3)dndp, 
n 

Cfe+i = Cfc - ^N k+ iAt 

An example of the time evolution of u(n, /3, t) and the corresponding con¬ 
sumption of the resource c{t) computed with the Algorithm 4 is shown in 
Fig- H Due to the mathematical equivalence of the two solutions, the states 
(snapshots) of u(n, /?, t) given by (22) are exactly equal to the states u(n, a , t) 
of the unmitigated exponential growth given by ([8]), provided a = (3. How¬ 
ever, the evolution of u(n , /3, t) proceeds at a different rate and freezes as 
c(t) —> 0. The accuracy of the Algorithm 4 can be eventually improved using 
the predictor-corrector technique. 


5 Populations with randomized migration 

Let a single metapopulation be nonuniformly distributed over some spatial 
domain at time t = 0. The whole habitat may be divided either naturally or 
virtually into a large number of elementary cells (sub-habitats) and a question 
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Fig. 3 Initial distribution uo(n, (3) (top-left), the final quasi-stationary distribution 
uo(n, fi, t m ax) (top-right), and the corresponding time-evolution of the resource c{t) (bot¬ 
tom). 


can be posed about the distribution of cell populations over their size at time t. 
The growth of these elementary cell populations could be described by one of 
the models from the preceding sections with growth parameters varying from 
cell to cell. If all individual members stay within their original cells or the 
whole metapopulation uniformly translates in space the phase-space methods 
developed earlier and the resulting distribution functions, obviously, apply 
without change. 

Allowing for the migration of individuals between the cells usually requires 
a PDE-based space-time dynamic law (e.g. Fischer equation). Such a law can¬ 
not be directly incorporated into the phase-space current, since the phase-space 
approach completely neglects the spatial ordering of the cells. Nevertheless, 
certain types of migration mechanisms can be described on the level of ODE’s 
and the corresponding phase-space continuity laws. 

Consider a stationary metapopulation, where the total number of individ¬ 
uals and the total number of sub-habitats (that play the role of populations 
in this case) do not change. Suppose that all individuals from time to time 
(with frequency /3) decide to emigrate from its native population so that each 
population is subject to the emigration rate /3n, proportional to its size. De¬ 
pending on the choice strategy of migrants as far as their new population is 
concerned one arrives at different ODE’s. 

For example, let each emigrant choose a new population (habitat) com¬ 
pletely at random. If viewed spatially, this is, obviously, a variant of super¬ 
diffusion with a random step size. Then, on a sufficiently large time scale, 
there will be the steady immigration rate 7 in all populations (since the total 
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emigration rate is constant). The emigration and immigration coefficients fi 
and 7 must be balanced to guarantee the steady state of the metapopulation: 


j t ^ ni = = °- 


(24) 


Let Y^i n i = N, then 7 = fin , where n = N/P , and P is the number of 
habitats. The dynamic law becomes 




i = 1,2,... P. 


(25) 


Notice that n = N/P is the total number of individuals divided by the number 
of habitats, which would be exactly the size of populations in all habitats if 


all individuals were evenly spread throughout. From (25) it follows that a 


population stops changing once it reaches this equilibrium size n. 

The distribution function u(n,t) is then found by solving the following 
Cauchy problem: 


du 

dt 


and is given by 


du 

)^n~fi u = Q, u(n, 0 ) 

= u 0 (n), 

(26) 

'P 1 uq (n— (n — n)e /3< ) , 

n <n, 


■ pt uo (n ), 

n = n, 

(27) 

•/ n u 0 (n + (n — rfije 13 *) , 

n>n. 



The middle branch of this solution shows that the number (density) of pop¬ 
ulations with n = n growth exponentially with time. As can be seen in Fig. [4] 
the upper and lower branches of the solution (27) represent the initial distri¬ 


bution that shrinks, respectively, from the left and from the right, towards n. 
This behaviour is indicative of a modifier of the Dirac delta function centred 
at n, which, indeed, is the limit of u(n,t) as t —> 00, since all populations will 
eventually have the same size n. An interesting feature of |27| ) is the preser¬ 
vation of the main features of the initial rt 0 (n) over time, albeit in a scaled 
(shrank) form, e.g., the double peaked shape shown in Fig. [4j 


6 Biased migration 

Consider populations with potentially unlimited growth and constant emigra¬ 
tion frequency. Let the choice of the new population by the migrants be biased 
towards populations with higher growth rates. The ODE describing this situ¬ 
ation is: 


dni 

dt 


otiUi - fini + 7 Oj, i = 1,... P, 


(28) 
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Fig. 4 Time evolution of the distribution function u(n,t) for populations with (uniform) 
migration given by the equation \27\ with 3 = 1.0 and n = 5.7. 


where on is the growth coefficient, /3 is the emigration frequency, and 7 is the 
immigration coefficient. Since migration does not change the total number of 
individuals (only the growth does), the following balance equation should be 
satisfied at all times: 


= 1^2 a u 


(29) 


reducing equation (28) to 
dn . 


Yu- 

— = (oi - P)m + cti/3^— 
at 2 ^ a j 


i,j = 1 


(30) 

(31) 

(32) 


The corresponding nonzero component of the phase-space current is 

J n — (a — (3)nu + a/3Ru , 

where 

R A ’ 

N = Jj nu(n,a,t) dadn, A = JJ au(n,a,t) dadn, 

and, since the a-component of the phase-space current is zero, A will stay 
constant in time. Thus, we arrive at the following weakly nonlinear continuity 
equation for the distribution u(n,a,t ): 

du du 

—- + (af3R + (a — /3)n)——I- (a — f$)u = 0, u(n, a, 0) = ^ 0(^7 a). (33) 

ot on 
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Fig. 5 Time evolution of the distribution function u{n , a, t) for populations with biased 
migration, see eq. ( |34| ). Top row - snapshots at t = 0.0,0.018, and 0.034; Bottom row - 
corresponding numerically integrated p(n, t ) showing the distribution of populations over 
their sizes. 


Assuming that R(t) is a given function of time we use the method of charac¬ 
teristics to obtain the implicit solution of this problem as 


u(n, a,t) = e ^ a ^uo (ne ^ — a/3£(a,t), , 

t 

o 


(34) 


However, R(t) and, hence, £(a,f) are functions of u as well, see eq. (32). Thus, 
we employ an iterative algorithm similar to the Algorithm 4. 


Algorithm 6 

— Given: £ 0 (a) = 0, uo(n, {3), N 0 , R 0 , and a sufficiently large 17. 

— For k = 0,1, 2,..., do: 


&+i(a) = &(<*) + Rke-^-^At, 

u k+1 (n,a) = e~^ a ~^ tk+1 UQ - a/3^ k+1 {a),ov) 


R k +1 = ^ JJ nu k+ i(n, a) dnda , 


A quick look at the implicit solution (34) reveals similarity with the case of 
exponential growth 0 - Indeed, as expected, the two solutions are identical for 
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(3 = 0, remain similar at low values of the emigration frequency, and become 
very different for larger values. One of the predictions of the metapopulation 
theory (Levin 1969; Eriksson et al. 2014) is that migration may prevent ex¬ 
tinction. In the present case, as follows from (30), none of the populations 
can disappear as long as there is at least one nonzero population somewhere. 
Yet, it does not mean that all populations are always growing in size. In the 
early stages some populations may decrease. From (30) we deduce that the 
condition on the growth of a population is 


^ f . 0^2 'N Hz 

v JJ Iv*? 


j3 > 0. 


(35) 


Obviously, this condition is satisfied for all populations with a* > f3. However, 
for a population with Oj < (3 it is easy to imagine a sufficiently large initial 
state rii(0) such that this condition is violated and the population size will 
decrease. Nevertheless, when such a decreasing rii becomes sufficiently small, 
the condition will become satisfied again and the population will resume its 
growth. 

An example of the evolution of distribution function for populations with 
biased migration computed with the Algorithm 6 is shown in Fig. [5} The initial 
distribution (the leftmost image and plot) is a Gaussian centred at n c = 20 and 
a c = 30. The emigration frequency is set at (3 = 50 in this simulation, which 
is higher than most of the growth coefficients a in the initial distribution. 
The first striking conclusion is that despite the immigration bias we do not 
observe an explosive growth of some particular (chosen) populations at the 
expense of others. Apparently, any such growth is efficiently mitigated by a 
higher emigration rate. 

The early-time transient decrease of some populations mentioned above 
(those with large initial n’s) causes the distribution function to contract slightly 
in the horizontal n-dimension, since the right side of the distribution moves to 
the left, while the left side moves to the right. This happens in the very early 
stages of the evolution (not shown). The same effect prevents the distribution 
from widening in the horizontal direction in the course of the evolution. In 
fact, as the Figure [5] shows, the distribution only keeps contracting. 

Thus, the main consequence of the large emigration frequency (3 for the 
late-time evolution appears to be the uniformity of size for populations with 
equal a’s and the uniformity of a’s for populations of the same size. Also, 
for large t’s the populations show a very strong correlation (almost direct 
proportionality) between the population size and its growth coefficient a (see 
the top-right image of Fig. [5]) . 


7 Growth rates depending on the distribution function 

The present approach is indispensable for dynamic models where the rate of 
growth of populations is a functional R(u, n) of both the size n of the local 
population and the distribution u of populations over the phase space. The 
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phase-space current vector will have as many nonzero components as there are 
time varying parameters in R. If all parameters, except n,, are time invariant, 
then the problem is described by the following equations: 


drii 

dt 

du 

at 


R(u,rii ), * = 1,2, ...P, 


0 

dn 


( uR(u , n )). 


(36) 

(37) 


Obviously, it is now impossible to solve the dynamic equations (361 without 
first solving the phase-space problem (371. 

Consider populations with unlimited growth discussed earlier. Suppose that 
one is able to influence the rate of growth of each population (e.g., by judicially 
watering and fertilizing each plant) proportionally to its current ‘weight’ in the 
phase space. Such a problem is described by the following equations: 


drii 

dt 


oti-\-Aot rii-\-An 


— Ql? 


i — Aot rin— An 


Xirii J J u(n,a,t) dnda ~ ain,iu(ni,ai,t)6, i = 1,2, ...P, 

(38) 

(39) 

where S = AAa An, i.e., we use the mid-point approximation for the integral 
in (38). The resulting approximate phase-space problem (39) is a variant of 
the Burgers equation (see e.g. Smoller 1994) whose characteristic equations 
are: 


du du 

— —I- 2aS nu— —I - aSu 2 = 0. u(n, a, 0) = uo(n, a), 
at on 


dt dn — du n 

— = 1, — = 2ao nu, — = —ad u , 

ds ds ds 


or 


dn = 2a6 nu dt, —a6 dt = 


du dn 
u 2 ’ 2 n 


du 


The solutions of the last two are: 


aSt -= C 2 , u\fn = C 3 , 


which allows rewriting and solving the first equation in (41) as 

OjTL 

—j= = 2a.dC?, dt, 2y/n(l — aS tu ) = C\. 

Jn 


Combining this result with the first equation in (42) we get 

1 


aSt- f(p)' 


(40) 


(41) 


(42) 


(43) 


u = 


p = 2y / n(l — aS tu). 


(44) 
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The initial condition requires that 


/(P)lt=0 = - 


u 0 (n,a)' 


(45) 


Hence, we choose 


fip) = - 


zto(p 2 /4, a) UQ{n(a5tu — ly ,a) : 


(46) 


and arrive at the following equation that implicitly defines the distribution 
u(n, a, t): 


or 


uo(n(aStu — l) 2 ,a) 
aS t uo(n(aS t u — l) 2 , a) + 1 


= (1 — aStu) Uo(n(l — a6t u) 2 , a). 


(47) 


(48) 


These representations are only valid in the smooth regime, i.e., prior to the 
development of shock. 

The fixed-point iterative algorithm that computes the distribution at a 
given time t may be formulated as follows (one can use any suitable norm for 
the computation of the current mismatch e k ): 

Algorithm 7 

— Given: uq( n,a), t, and e. 

— For k = 0,1, 2,, while e k > e, do: 

Uo(n(aS tuh(n, a) — l) 2 , a) 


k+1 ^ ’ ^ aSt uo(n(aS tuh(na) — l) 2 , a) + 1 
e k = \\u k+ i(n,a) - u k (n,a)\\ . 

Figure [6] illustrates the application of this algorithm (with e = 10~ 4 and 
5 = 0.1). The initial distribution is the same as in the example of the previous 
section. As expected from the non-viscous Burgers equation the solution devel¬ 
ops a shock (notice the almost vertical front in the plot at the bottom-right of 
Fig# Even though we have used the distribution computed at the previous 
time instant as the initial guess for the next time instant, the convergence of 
the Algorithm 7 slows down as time grows. Around t = 0.1, as one approachers 
the shock in this example, the convergence of the Algorithm 7 breaks down 
completely: initially the error e k decreases, however, it never reaches the tol¬ 
erance of e = 1CT 4 and even begins to increase as iterations continue. 

Obviously the development of shock and the related breakdown of the 
Algorithm 7 does not mean that the populations stop growing. Rather, it is 


a result of the midpoint approximation of the integral in the equation (381. 


Keeping that in mind, we see that at least in the early stages of evolution 
there will emerge a sharp upper bound on the sizes of populations and many 
populations with the average growth factor will have that size. 
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Fig. 6 Time evolution of the distribution function u[n,a,t) for populations with growth 
rates depending on u, see eq. ( |47[ ). Top row — snapshots at t = 0.0, 0.054, and 0.09; Bottom 
row — corresponding numerically integrated p(n, t) showing the distribution of populations 
over their sizes. 


8 Conclusions and possible extensions of the model 

The phase-space approach to multi-population dynamics presented above al¬ 
lows computing the time evolution of the distribution function that approxi¬ 
mates the multi-dimensional histogram of populations. Although we have con¬ 
sidered some trivial stochastic processes in the sections about migration, the 
main value of this approach is in the fast and easy analysis of causal de¬ 
terministic trends in biological data. For example, in the case of precision 
agriculture (Hautala and Hakojarvi 2011), this method could be employed 
to predict and, perhaps, control the distribution of plant sizes at the time 
of harvest. To this end, instead of introducing stochastic variations due to 
the environment, models considered above can be easily extended to include 
deterministic (measured) time-varying growth coefficients and environmental 
factors. Another straightforward extension could be the inclusion of the source 
term in the conservation law ([5]) reflecting the emergence of new populations. 

The case of populations competing for limited resources could serve as a 
mathematical model for recovering the distribution of populations from the 
resource consumption curves. For example, in the analysis of seed quality 
it is common to monitor the oxygen consumption during the germination 
process (Bewley et al. 2013, van Duin and Koenig 2001, van Asbrouck and 
Taridno 2009), with the main consumers of oxygen being the growing popula¬ 
tions of mitochondria present in every seed. Our method provides the neces- 
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sary link between the distribution of the seed parameters, including the initial 
number of active mitochondria, and the measured oxygen uptake curves for a 
batch of seeds germinating in a closed container. 

It is easy to imagine that the phase-space dynamics would be even more 
exciting if the phase-space current had several nonzero components, e.g., with 
the time varying growth/environmental factors or with the classical predator- 
prey problem. The corresponding continuity equations, however, would have 
to be solved numerically (Leveque 2002). 

Although we were able to analyze two special migration problems, the 
general case of linear coupling between ODE’s so far appears to be intractable 
with our methodology. This has mainly to do with the nonlinearity of the 
mapping between a given vector and its histogram or distribution. Indeed, 
the histogram of the sum of two vectors is not the sum of the histograms of 
these vectors. Instead, we should, probably, focus on the phase-space coupling 
considered in the last section. 


Acknowledgements and ethical statement 

The authors are grateful to Prof. A. van Duijn (Leiden University and Fytago- 
ras BV) and Dr. S. Hille (Leiden University) for introduction into the field of 
plant physiology and the problem of seed germination, and to Dr. H. van Doom 
(HZPC Holland BV) for introduction to the problems of precision agriculture. 
The authors declare no conflicts of interest. 


References 

1. Begon M, Townsend CR, and Harper JL (2006) Ecology: From Individuals to Ecosys¬ 
tems. Blackwell Publishing. 

2. Bewley JD, Bradford KJ, Hilhorst HWM, and Nonogaki H (2013) Seeds: Physiology of 
development, germination and dormancy. 3rd Edition, Springer, New York, Heidelberg, 
Dordrecht, London. 

3. Carrie C, Murcha MW, Giraud E, Ng S, Zhang MF, Narsai R, and Whelan J (2012) 
How do plants make mitochondria? Planta: An International Journal of Plant Biology, 
10.1007/ s00425-012-1762-3. 

4. Collet JF and Goudon T (2000) On solutions of the Lifshitz-Slyozov model. Nonlinearity 
13:1239-1262. 

5. Edelstein-Keshet L (2005) Mathematical Models in Biology. SIAM. 

6 . Eriksson A, Elias-Wolff F, Mehlig B, and Manica A (2014) The emergence of the rescue 
effect from explicit within- and between-patch dynamics in a metapopulation. Proc. R. 
Soc. B 281:20133127. 

7. Fisher RA (1930) The genetical theory of natural selection. Oxford University Press. 

8 . Hautala M and Hakojarvi M (2011) An analytical C3-crop growth model for precision 
farming. Precision Agric. 12:266—279. 

9. Kampman R and Wagner R (1991) Materials Science and Technology - A Comprehensive 
Treatment. Vol. 5, VCH, Weinheim, Germany. 

10. Lame ryot P (2002) The Lifshitz-Slyozov-Wagner equation with conserved total volume. 
SIAM J. Math. Anal. 34(2):257-272. 

11. Leveque RJ (2002) Finite volume methods for hyperbolic problems. Cambridge Univer¬ 
sity Press, Cambridge. 



18 


Neil V. Budko, Fred J. Vermolen 


12. Levin R (1969) Some demographic and genetic consequences of environmental hetero¬ 
geneity for biological control. Bull. Entomol. Soc. Am. 15:237-240. 

13. Lifshitz IM and Slyozov VV (1961) The kinetics of precipitation from supersaturated 
solid solutions. J. Phys. Chem. Solids, 19:35-50. 

14. Lin Y, Berger U, Grimm V, Huth F, and Weiner J (2013) Plant Interactions Alter the 
Predictions of Metabolic Scaling Theory. PLoS ONE 8(2):e57612. 

15. Merchant DJ, Kuchler RB, and Munyo W (1960) Population dynamics in suspension 
cultures of an animal cell strain. Journal of Biochemical and Microbiological Technology 
and Engineering, 11:253-266. 

16. Myhr OR and Grong 0 (2000) Modelling of non-isothermal transformations in alloys 
containing a particle distribution. Acta Materialia, Vol. 48(7): 1605-1615. 

17. Newman K, Buckland ST, Morgan B, King R, Borchers DL, Cole D, Besbeas P, 
Gimenez O, and Thomas L (2014) Modelling Population Dynamics: Model Formula¬ 
tion, Fitting and Assessment using State-Space Methods. Series: Methods in Statistical 
Ecology, Springer Science+Business Media, New York. 

18. Shonkwiler RW and Herod J (2009) Mathematical Biology: An introduction with Maple 
and Matlab. 2nd Edition, Springer Dordrecht, Heidelberg, London, New York. 

19. Smoller J (1994) Shock waves and reaction-diffusion equations. Springer-Verlag, New 
York, Second Edition. 

20. van Asbrouck J and Taridno P (2009) Using the single seed oxygen consumption mea¬ 
surement as a method of determination of different seed quality parameters for com¬ 
mercial tomato seed samples. Asian Journal of Food and Agro-Industry:S88—S95. 

21. van Duijn A and Koenig W (2001) Measuring metabolic rate changes. Patent EP1134583 
(Al) WO0169243 (Al) US2004033575 (Al) CA2403253 (Al) DE60108480T. 

22. Velazquez J, Garrahan JP, and Eichhorn MP (2014) Spatial Complementarity and the 
Coexistence of Species. PLoS ONE 9(12):ell4979. 

23. Wagner C (1961) Theorie der Alterung von Niederschlagen durch Umlosen (Ostwald- 
Reifung). Z. Elektrochem., 65:581-591. 

24. Xu L, Zhang F, Zhang K, Wang E, and Wang J (2014) The Potential and Flux Land¬ 
scape Theory of Ecology. PLoS ONE 9(l):e86746. 



